Rcwl_MC3

Introduction

The Multi-Center Mutation Calling in Multiple Cancers project (MC3) pipeline was developed by TCGA to generate a comprehensive encyclopedia of somatic mutation calls by applying an ensemble of seven mutation-calling algorithms with scoring and artifact filtering. More details can be found in this paper: Scalable Open Science Approach for Mutation Calling of Tumor Exomes Using Multiple Genomic Pipelines

The mc3 pipeline is available at https://github.com/OpenGenomics/mc3. All required software have been deployed in cloud with docker.

The pipeline has been imported and contained in the RcwlPipelines pacakge, which contains two major steps (markID step was removed):

• Call variants by 7 pipelines
• Merge VCF and convert to MAF

Here is the short summary.

suppressPackageStartupMessages(library(Rcwl))
library(RcwlPipelines)

data(mc3)
short(mc3)

inputs:
- tumorID
- normalID
- tumor
- normal
- bed_file
- centromere
- cosmic
- dbsnp
- refFasta
- vepData
outputs:
- outmaf
- outvcf
steps:
- call_variants
- convert

plotCWL(mc3)

1. call_variants

callVar <- readCWL(runs(mc3)$call_variants)
plotCWL(callVar)

2. Convert

conv <- readCWL(runs(mc3)$convert)
plotCWL(conv)

Prepare data

Testing somatic mutation data can be download from: https://github.com/genome/somatic-snv-test-data.

1. Input list inputList. The tumorID/normalID must be consistent with SM from BAM read group.

inputList <- list(tumorID=list(test="NA12892"),
                  normalID=list(test="NA12878"),
                  tumor=list(test="data/tumor.bam"),
                  normal=list(test="data/normal.bam"))

2. Parameter list paramList.

paramList <- list(bed_file="/mnt/lustre/users/qhu/software/mc3/data/gaf_20111020+broad_wex_1.1_hg19.bed",
                  centromere="/mnt/lustre/users/qhu/software/mc3/data/centromere_hg19.bed",
                  cosmic="/mnt/lustre/users/qhu/software/mc3/data/hg19_cosmic_v54_120711.vcf.gz",
                  dbsnp="/mnt/lustre/users/qhu/software/mc3/data/dbsnp_134_b37.leftAligned.vcf.gz",
                  refFasta="/rpcc/bioinformatics/reference/current/human_g1k_v37.fa.gz",
                  vepData="/home/qhu/.vep/")

Run MC3 pipeline

res <- runCWLBatch(mc3, outdir = "output/mc3",
                   inputList = inputList, paramList = paramList,
                   BPPARAM = BatchtoolsParam(workers = 1, cluster = "sge",
                                             template = "/rpcc/bioinformatics/sge_centos7.tmpl",
                                             resources = list(threads = 2,
                                                              queue = "centos7.q")))

The final VCF was filtered and merged from the outputs of different pipelines and annotated by VEP. The converted MAF file was also generated.

dir("output/mc3/test")

[1] "merged.vep.vcf" "vep.maf"       

vcf <- read.table("output/mc3/test/merged.vep.vcf", sep="\t")
head(vcf)

  V1       V2          V3 V4 V5  V6   V7
1 21 10400299           .  A  T 0.0 PASS
2 21 10400380           .  C  T   . PASS
3 21 10402435   rs2948877  G  A   . PASS
4 21 10402715           .  G  A 0.0 PASS
5 21 10402795 rs148043841  G  T   . PASS
6 21 10402985           .  G GA   . PASS
                                                                                                                                                                                               V8
1                                               CENTERS=RADIA|VARSCANS|MUSE|SOMATICSNIPER;CSQ=T|intergenic_variant|MODIFIER|||||||||||||||rs370695467|1||||1|SNV|1|||||||||||||||||||||||||||||||
2                                                          CENTERS=SOMATICSNIPER|RADIA|VARSCANS|MUSE;CSQ=T|intergenic_variant|MODIFIER||||||||||||||||1||||1|SNV|1|||||||||||||||||||||||||||||||
3 CENTERS=MUSE|RADIA|VARSCANS|SOMATICSNIPER;CSQ=A|intergenic_variant|MODIFIER|||||||||||||||rs2948877|1||||1|SNV|1||||||||||||||||A:0.3626|A:0.4266|A:0.3228|A:0.3075|A:0.2913|A:0.4346||||||||||
4                                                               CENTERS=RADIA|VARSCANS|MUSE;CSQ=A|intergenic_variant|MODIFIER|||||||||||||||rs2948878|1||||1|SNV|1|||||||||||||||||||||||||||||||
5                                                             CENTERS=MUSE|RADIA|VARSCANS;CSQ=T|intergenic_variant|MODIFIER|||||||||||||||rs373568457|1||||1|SNV|1|||||||||||||||||||||||||||||||
6                                                          CENTERS=VARSCANI*|PINDEL;CSQ=A|intergenic_variant|MODIFIER|||||||||||||||rs375209288|1||||1|insertion|1|||||||||||||||||||||||||||||||
        V9           V10            V11
1 GT:DP:AD 0/0:140:140,0   0/1:92:71,20
2 GT:DP:AD 0/0:160:160,0  0/1:117:99,18
3 GT:DP:AD 0/0:167:167,0  0/1:124:97,27
4 GT:DP:AD 0/0:145:145,0  0/1:117:97,20
5 GT:DP:AD 0/0:163:161,2 0/1:127:108,19
6 GT:DP:AD   0/0:88:88,0    0/1:82:75,7

Session information

sessionInfo()

R version 3.5.2 Patched (2018-12-31 r75935)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS release 6.4 (Final)

Matrix products: default
BLAS: /home/qhu/usr/R-3.5/lib64/R/lib/libRblas.so
LAPACK: /home/qhu/usr/R-3.5/lib64/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] RcwlPipelines_0.99.0 testthat_2.0.1       BiocStyle_2.10.0    
 [4] rmarkdown_1.11       knitr_1.21           usethis_1.4.0       
 [7] devtools_2.0.1       Rcwl_0.99.14         S4Vectors_0.20.1    
[10] BiocGenerics_0.28.0  yaml_2.2.0          

loaded via a namespace (and not attached):
 [1] colorspace_1.4-0    rprojroot_1.3-2     fs_1.2.6           
 [4] rstudioapi_0.9.0    roxygen2_6.1.1      remotes_2.0.2      
 [7] xml2_1.2.0          pkgload_1.0.2       jsonlite_1.6       
[10] workflowr_1.2.0     shiny_1.2.0         DiagrammeR_1.0.0   
[13] BiocManager_1.30.4  readr_1.3.1         compiler_3.5.2     
[16] backports_1.1.3     assertthat_0.2.0    lazyeval_0.2.1     
[19] cli_1.0.1           later_0.8.0         visNetwork_2.0.5   
[22] htmltools_0.3.6     prettyunits_1.0.2   tools_3.5.2        
[25] igraph_1.2.4        gtable_0.2.0        glue_1.3.0         
[28] dplyr_0.8.0.1       batchtools_0.9.11   rappdirs_0.3.1     
[31] Rcpp_1.0.0          rgexf_0.15.3        xopen_1.0.0        
[34] debugme_1.1.0       xfun_0.5            stringr_1.4.0      
[37] ps_1.3.0            mime_0.6            XML_3.98-1.17      
[40] scales_1.0.0        hms_0.4.2           promises_1.0.1     
[43] RColorBrewer_1.1-2  memoise_1.1.0       gridExtra_2.3      
[46] ggplot2_3.1.0       downloader_0.4      rcmdcheck_1.3.2    
[49] stringi_1.3.1       highr_0.7           Rook_1.1-1         
[52] desc_1.2.0          checkmate_1.9.1     pkgbuild_1.0.2     
[55] BiocParallel_1.16.6 rlang_0.3.1         pkgconfig_2.0.2    
[58] commonmark_1.7      evaluate_0.13       purrr_0.3.0        
[61] htmlwidgets_1.3     processx_3.2.1      tidyselect_0.2.5   
[64] plyr_1.8.4          magrittr_1.5        bookdown_0.9       
[67] R6_2.4.0            base64url_1.4       pillar_1.3.1       
[70] withr_2.1.2         tibble_2.0.1        crayon_1.3.4       
[73] viridis_0.5.1       progress_1.2.0      grid_3.5.2         
[76] data.table_1.12.0   git2r_0.24.0        callr_3.1.1        
[79] influenceR_0.1.0    digest_0.6.18       xtable_1.8-3       
[82] tidyr_0.8.2         httpuv_1.4.5.1      brew_1.0-6         
[85] munsell_0.5.0       viridisLite_0.3.0   sessioninfo_1.1.1